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Abstract 

We use a boolean cellular automaton model to describe the diffusion limited dynamics 
of the irreversible reaction A + A^A + S on a, ID lattice. We derive a set of equations 
for the dynamics of the empty interval probabilities from which explicit expressions for 
the particle concentration and the two point correlation can be obtained. It is shown 
that the long time dynamics is in agreement with the off-lattice solution. The early time 
behavior, however, predicts a slower decay of the concentration. 
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I. INTRODUCTION 



In a recent paper ben-Avraham et al. treated the ID dynamics of the diffusion 
limited reaction A + A — > A + S. The analysis led to an infinite hierarchy of differential 
equations (HDE) for the time evolution of the set of quantities {Ek(t)}, where Ek(t) 
denotes the probability of finding an interval of k contiguous empty lattice sites (c.f. eq. 
(2.2) and (2.5) in |l| ). Ben-Avraham et al. were primarily interested in the continuum 
limit for the probability densities Ek(t)/(Ax) k as the lattice spacing Ax is allowed to 
vanish. Here we rederive the HDE taking as a starting point a detailed cellular automaton 
description using Boolean occupation numbers Instead of going to the continuum 

limit, we solve the HDE for a discrete lattice and obtain explicit expressions for the 
concentration c(t) = (1 — Ei(t)) /Ax and the two-point nearest neighbour correlation. 
We consider both the case of an initially fully occupied lattice and the case in which the 
initial probabilities Ek(0) are given by the random "geometrical" distribution (1 — p) k , 
where p < 1 is the initial probability of a single site being occupied. Some exact solutions 
to related problems and a special case can be found in ||. 

In section [TT], we derive the HDE starting from the cellular automaton dynamical rule. 
The next section is a brief reminder of the off-lattice solution for the HDE obtained in 
JT|. In section [TV], we solve the solution of the HDE for a discrete lattice and discuss 
the transient dynamics of the concentration. As expected, the long time behavior is the 
same as that found in ffl. However, the early time behavior turns out to be significantly 
different. 

II. DERIVATION OF THE HDE 

Consider a ID lattice with N sites. Each lattice site i is characterized by a boolean 
occupation number n, = 1 or depending on whether it is occupied by a single particle 
(A) or empty (S). At each time step 5t, we choose randomly a site i in the lattice via 
the boolean stochastic parameter . This parameter is equal to one for the chosen site 
and for all other sites. Simultaneously, a second boolean variable is used to select the 
left (£l = 1) or the right neighbour site (£x = 0) with equal probability 1/2. If the site 
i is occupied, the particle will hop to the chosen neighbour site with a rate hp given by 
the mean value of the stochastic boolean parameter £q. In this case, if the neighbour 
site is empty, the particle will occupy it vacating the site i. We express this by including 
the loss terms 

-^{V^U^nMil-n^it)) (1) 
(1 - &(*)ni(t)(l - n l+1 (t)) (2) 

in the dynamical rule for rij(t). If the neighbour site is filled, the particle will "react" 
with it and instantaneously disappear. Again, the occupation number at site i will be 
decreased from 1 to 0. Thus, we have the reactive loss terms: 
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-£#(*KL(*)fr(*M*)ni-i(f) (3) 
-$(t) (l-^ L (t))^ D (t) ni (t)n i+1 (t) (4) 

Finally, as an empty site i can only be occupied by hopping from a particle at a neighbour 
site, one has the two gain terms: 

+d +1) (t)Ut)Ut)n i+1 (t)(l-n i (t)) (5) 

+dr 1} w (i - £D(t)"i-i(*)(i - (6) 

Clearly, the dynamics described above corresponds to the particular implementation of 
the reaction A + A — > A + S given in |L|. Adding up all the contributions ([[])- (||), we 
obtain the following dynamical rule: 

m(t + 5t) = m(t) - $(t) £D(t)ni(t) 

+dT 1} W (1 - &(*)n»_i(t) (1 - n, t (t)), i = 1, . . . , N. (7) 

In the first and last equation (0), we set n (t) = n N+1 (t) = 0. We can rewrite (0) by 
using the complementary occupation numbers Sj(t) = 1 — nj(i): 

+ = 8*® + ^(t) ut) - izVHt) + &\m D (t>i(t) 

+^- 1 \t)(l-^ L (t))Ut)s l -i(t) Sl (t) 

(t) tD(t) Si (t) s i+1 (t), i = 1, . . . , N. (8) 



This is a more convenient form, since, as it turns out, the evolution law for a string 

Tl+fc-l „ 

Lj=i 6 J 



of k consecutive empty sites YViLi 1 s j involves only products of contiguous occupation 



numbers 

^_|_^- ^ j-j-fo 1 i-\-k 1 

n Sj (t+At)= n s j (t)+^\t)ut)ut) n 

j=i j=i j=i+l 

i+k-2 

j=i 

-fccwwoTElUw +e{^ ) (*)a(*)fo(on«i(*) (9) 

Let us take the time step St equal to 1/N, implying that each site has been visited once 
on average after one time unit. If we now average (Q) over an ensemble of realizations 
for a given initial configuration, we obtain in the thermodynamic limit iV — > oo: 
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d -§ = y (^i + E ^ - 4 ^ + E ^ + E Ui) , (io) 



where E\{t) = Ujti Sj(t). For a single site (k = 1) the corresponding equation 

^ -k D (l- 1E\ + i^- 1 + \e\) (11) 



di " V 1 2 z 2 

is obtained by averaging the dynamical rule (H). If we perform the average over re- 
alizations and translationally invariant initial conditions, we get the following hier- 
archy of differential difference equations for the evolution of the averaged products 

E k (t) = (m- 1 s 3 m 



dE k 



dt 



k D {E k+1 - 2E k + E k - X ) , fc = 1,2,... (12) 



with the boundary condition E (t) = 1. As described in [[I]], the rhs of QT^ ) represents 
the net flux due to particle diffusion into and out of an empty site interval, whereas the 
effect of reaction enters through the boundary condition. 



III. OFF-LATTICE SOLUTION 

Following ref. Q, we set the hopping rate k D /2 to either of both sides equal to 
D/(Ax) 2 , where Ax is the lattice spacing. On long length and time scales this yields 
normal diffusion with a diffusion coefficient D. With this definition, eqs. (0) are 
identical to those derived in 0. If we now let Ax — > 0, eqs. (O) become 



dE(x,t) n „ d 2 E(x,t) 

dt ~ lD dx 2 ■ [L6) 

with the boundary conditions £7(0, t) = 1 and E(oo,t) = 0. In this limit, the concentra- 
tion (number of particles per unit length) is expressed as 



c{t) = _dE( x , t ) 



dx 



(14) 



x=0 



Thus, one can determine the time dependence of the concentration by solving (|13|) with 
the boundary conditions given above (see [|TJ for details). For the special case of an 
initially random particle distribution with a concentration cq, one has 



c(t) [8c 2 Dt\^ 



C \ TT 

and 



+ o(c 2 Dt) as t -> (15) 



c{t) — > r as t — > oo (16) 

(27rDt)5 
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IV. SOLUTION FOR A DISCRETE LATTICE 



We now proceed to solve the HDE (|I2|) for the case in which one has an initially full 
lattice, i.e. E k (0) = for all k > 1. To begin with, we absorb the rate constant kn into 
the time scale by introducing the dimensionless time variable t = hpt. The hierarchy 
then reads: 

^■ = E k+1 -2E k + E k . u k = 1,2,... (17) 

GST 

Next we apply the Laplace transform to both sides of ([II]) and obtain the homogeneous 
difference equation 

E k+1 -(2 + s)E k + E k ^ = 0, k>l, (18) 
where E k (s) = L T ^ S {E k (r)} = / exp(— st) E k (r) dr. The boundary condition is 

_ JO 

given by E (s) — 1/s. This second-order difference equation is solved with the ansatz 
Ek(s) = X k (s). This leads to the quadratic equation 

A 2 -(2 + s)A + l = (19) 

which has the two solutions 

A ± = i+l^HE (20) 

The general solution of (|1|) is obtained as a linear superposition of and A fc : 

E k (s) = A(s) X k _(s) + B(s) A+(s) (21) 



For k — * oo, the physically acceptable solution of (|T8|) must satisfy the implicit boundary 
condition E^s) = 0. To avoid the divergence of the second term in (|2l|) , we must 
therefore set B(s) = 0. Using the other boundary condition at k — 0, we find A(s) = 1/s. 
Thus, we have: 



BtW = i ( £±lz^±l ) (22) 

Clearly, the most interesting quantity is Ei(s), whose inverse Laplace transform Ei(r) 
is the probability that a randomly chosen site be empty. By virtue of a theorem 0, the 
inverse transform L~}^ T {E k (s)} is given by the integral 

v k (r')dr', k = 1,2, (23) 

where 
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Mt) = L-X { | i±2^EIE\\ = t exp(-2r) 4(2.) (M) 



(see e.g. 0, page 379]), where 



I n (x) = > — — — r 25 



are the modified Bessel functions. In particular, one has 

exp(-2r')/i(2r') _, , 



jo r 



A. Asymptotics for early times 

For sufficiently short times (r <C 1), we can use (^Sj) and the series expansion of 
the exponential function exp (— 2r') to expand the integrand in (|26|) in powers of r'. 
Neglecting terms of order o(r' 3 ), performing the integration and undoing the time scaling, 
we obtain: 

E 1 (t) = k D t-k 2 D t 2 + ^k 3 D t 3 + o(t A ) (27) 

o 

The particle concentration is 

c(t) = ^= =c \l- 2<*Dt + KDH 2 + o{clDH 3 )] , (28) 

Ax Ax L J 

where Co = 1/Ax. This is in clear disagreement with the decay law fll5"l). 

For sufficiently short times, we can neglect the effect of large clusters of vacant sites, 
since the chain is initially full. This is done by setting Ek = for k larger than a certain 
cutoff size k c in the truncation hierarchy flTT|). If we set k c = 1, we obtain the differential 
equation 

-^ = l-2^(r). (29) 

(XT 

The solution reads 

E 1 (t) = \ (1 - exp (-2r)) = r - r 2 + ^r 3 + o(r 4 ) (30) 

Setting t = k D t, we see that the early times expansion (^7j) is reproduced correctly up to 
the quadratic term. The discrepancy between ( p7|) and the off-lattice solution arises due 
to the finite propagation velocity of a local perturbation in concentration, as opposed 
to the infinite propagation velocity characteristic of diffusion. Thus, we expect a slower 
decay of the concentration c(t) on the lattice. 
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B. Long time asymptotics 



For large times (r ^> 1), we write Ei(t) as follows: 



, exp (-2r')/i(2r') , , 
E 1 (r) = A- — ^ ^ J -dr'. (31) 



dr', (32) 



where A is the definite integral 

exp(-2r , )/i(2r') 



r' 



which is equal to 1 |[ page 236]. The integrand in the second term of ( |3~TD can be 
expanded using the asymptotic form 

«•> -3^-5 -(?))■ < 33 ' 

for large x ( see e.g. ]7|, page 489] ). Thus, we get 

1 r°° a 3 r°° 5 5 

= 1 -^= + —^ + 0(^-4). (34) 
This is in agreement with the long time asymptotics of the continuum limit solution fllBD . 

C. Two-point correlation 

The explicit expression for the two interval probability 

r exp (-2r')I 2 (2r') , , 
E 2 (r) = 2 — ^ — J -dr'. (35) 



o 



can be used as a starting point to compute the asymptotics of the two-point nearest 
neighbour correlation c^ 2 '(t) = P 2 /(Ax) 2 , where P 2 = 1 — 2E\ + E 2 is the probability 
of finding two contiguous particles in the lattice. For early times one gets c^ 2 \t) = 
c 2 [1 - c 2 Dt + o(c*DH% whereas for long times c^(t) « l/c VS2^(Dt) 3 / 2 . 



D. Solution for an arbitrary initial concentration 

In this case, we average over all possible initial configurations of the lattice with a 
given concentration cq = p/Ax (p < 1). The initial conditions for the empty interval 
probabilities now read 

E k (0) = (l-p) fc , fc = l,2,... (36) 
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The boundary conditions are the same as in the case of an initially full lattice. If we 
now apply the Laplace transform to fll7), we obtain 

E k+1 - (2 + s)E k + E k ^ = -(1 - p)\ k > 1, (37) 

which differs from ([Tj|) by the inhomogeneity on the rhs. Like in the theory of ordinary 
differential equations, the general solution of ([57]) can be written as the sum of the general 
solution (PH) for (|T8|) and a particular solution which we seek in the form 

E P k a \s) = {l-p) k C{s). (38) 

Inserting this ansatz in (p7|), we find 

C(s) 



(1-P) 



(1 - p)s- p 2 



(39) 



The boundary condition for k — > oo again imposes B(s) = 0. From the boundary 
condition for E (s) we obtain 



A(s) 



1 



1-P 



s (1 — p)s — p 2 



Putting (|3^) and fl40|) into the equation 

£ fc ( s ) = A( s )A"( s ) + £T 

we find 

1-p 



1 



s (l-p)s-p 2 



s + 2 - y/s 2 + 4s 



+ 



fc+1 



;i - p)s - p 2 



k = 0,1, 



(40) 



(41) 



(42) 



We can now use the convolution theorem for the Laplace transform to invert ([42]). 
This yields 



EJt) = k 



+ 



exp (-204(20^ 



t' 



[1-pf-k 



exp ^— 



2 + ^ 



r')/ fe (2r') 



dr 7 



exp 



(43) 



However, we can study the asymptotics directly from ( f4"2"D by making use of the Tauberian 
theorems ||, which allow us to determine the behaviour of Ei{t) for r — ► and r — > oo 
by inverting respectively the series expansion of Ei(s) around s = oo and s = 0. For 
small s one has 



1 1 1 1 -P 

S y/S 2 p 



(44) 
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from which we get 



E l « 1 



for t — > oo. 



7TT 



(45) 



Again, this is in agreement with (|IJ). The long time asymptotics does not depend on 
the initial concentration cq, suggesting that a universal behavior also takes place on a 
finite lattice. In the opposite limit we have 



s s s a 
Inverting term by term the rhs of fl4"BD, we get 

E 1 = 1 - p + p 2 r - p2(1 + p) r 2 + o(r 3 ) 



(46) 



(47) 



leading to 



c(t) = c 



- i/,J(/M 1! r ( ;/; 2 / 2 • ^r;;/; 1 / 1 ; 



1 - -c^t + 2 

p p 5 



(4J 



Thus, the discrepancy with (|TJ) appears to be robust. We can again compare the exact 
early time expansion (|47[) with the solution of the differential equation fl29 ) with the 
initial condition -E'i(O) = 1 — p, which is given by 



- + - - p exp (-2r) = 1 - p + {2p - l)r + (1 - 2p)r 2 + o(r 3 ). 



(49) 



As expected, the approximation based on the neglect of empty intervals becomes worse 
as p decreases. 
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